Current Year Sea Surface Temperature Anomaly Maps

Aerial views of global temperature patterns

Author
Affiliation

Gulf of Maine Research Institute

Published

August 29, 2022

Code
# Packages
library(lubridate) # date support
library(raster)  # raster support
library(rnaturalearth) # coastline polygons
library(sf) # simple feature support
library(stars) # raster plotting with ggplot
library(gmRi) # styling for GMRI
library(scales) # axis labels
library(here) # project navigation
library(janitor) # data cleaning
library(patchwork) # multiplot arrangement
library(tidyverse) # data wrangling and plotting
library(xaringanExtra) # tab panels
library(ggthemes) # theme support



# Support Functions
suppressPackageStartupMessages(source(here("R/oisst_support_funs.R")))
suppressPackageStartupMessages(source(here("R/temp_report_support.R")))
path_fun <- os_fun_switch(mac_os = "mojave")


# File Paths
mills_path <- path_fun("mills")
oisst_path <- path_fun("res", "OISST/oisst_mainstays")


# Set ggplot theme for figures
theme_set(theme_bw() + 
  theme(
    # Titles
    plot.title = element_text(hjust = 0, face = "bold", size = 14),
    plot.subtitle = element_text(size = 9),
    plot.caption = element_text(size = 7.2, margin = margin(t = 20), color = "gray40"),
    legend.title = element_text(size = 9),
    legend.text = element_text(size = 7.5),
    # Axes
    axis.line.y = element_line(),
    axis.ticks.y = element_line(), 
    axis.line.x = element_line(),
    axis.ticks.x = element_line(), 
    axis.text = element_text(size = 11),
    # Facets
    strip.text = element_text(color = "white", 
                              face = "bold",
                              size = 11),
    strip.background = element_rect(
      color = "white", 
      fill = "#00736D", 
      size = 1, 
      linetype="solid")
    )
  )

# Polygons for mapping
new_england <- ne_states("united states of america", returnclass = "sf")
canada      <- ne_states("canada", returnclass = "sf")
world       <- ne_countries(returnclass = "sf")
greenland   <- ne_states(country = "greenland", returnclass = "sf")

# Adding Logo to plot
logo_path <- paste0(system.file("stylesheets", package = "gmRi"), "/gmri_logo.png")
lab_logo <- magick::image_read(logo_path)
Code
# Set parameters:
params <- list(
  region = "apershing_gulf_of_maine",
  in_C = TRUE
)


# File paths for various extents based on params$region
region_paths <- get_timeseries_paths(region_group = "gmri_sst_focal_areas", 
                                     mac_os = "mojave")
  
# clean up the name
tidy_name <- str_replace_all(params$region, "_", " ") %>% str_to_title() %>% str_replace(pattern = "Of", replacement = "of")
tidy_name <- str_remove(tidy_name, "Aak ")
tidy_name <- str_remove(tidy_name, "Apershing ")
tidy_name <- str_remove(tidy_name, "Cpr ")


# Current Year
current_yr <- year(Sys.Date())

# Use panelset
use_panelset()
Code
# Turn on the stylesheet:
gmRi::use_gmri_style_rmd(css = "gmri_rmarkdown.css")
Code
# Turn on the gmri font for plots
showtext::showtext_auto()

1: Global Temperature Maps

The 2022 global sea surface temperature anomalies have been loaded and displayed below to visualize how different areas of the ocean experience swings in temperature.

For perspective on where excess heat is distributed around the world, here first are maps of anomalies at a global perspective:

Code
# Access information to netcdf on box
nc_year <- "2022"
anom_path <- str_c(oisst_path, "annual_anomalies/1982to2011_climatology/daily_anoms_", nc_year, ".nc")

# Load 2022 as stack
anoms_2022 <- stack(anom_path)

Annual Average

Code
# Make Annual Average
ann_anom_ras <- calc(anoms_2022, fun = mean, na.rm = T)

# Color scale title
sst_lab <- "Sea Surface Temperature Anomaly"

# Color limit for palettes
temp_limits <- c(-2, 2)
temp_breaks <- c(temp_limits[1], temp_limits[1]/2,  0, temp_limits[2]/2, temp_limits[2])
temp_labels <- str_c(c(str_c("< ", temp_limits[1]), temp_limits[1]/2, 0, temp_limits[2]/2, str_c("> ", temp_limits[2])), "\u00b0C")

# Make st object
ann_anom_st <- st_as_stars(rotate(ann_anom_ras)) 

# Build Map
ggplot() +
  geom_stars(data = ann_anom_st) +
  geom_sf(data = world, fill = "gray30", color = "white", size = .25) +
  scale_fill_distiller(
    palette = "RdBu", 
    na.value = "transparent", 
    limit = temp_limits, 
    oob = scales::squish,
    breaks = temp_breaks, 
    labels = temp_labels) +
  guides("fill" = guide_colorbar(
    title = sst_lab, 
    title.position = "right", 
    title.hjust = 0.5,
    barheight = unit(3, "inches"), 
    frame.colour = "black", 
    ticks.colour = "black")) +  
  coord_sf(expand = F) +
  map_theme() +
  theme(legend.position = "right", 
        legend.title = element_text(angle = 90)) +
  labs(title = str_c("Global Temperature Anomalies - ", nc_year, " through: ", Sys.Date()))

Winter

Code
# Get previous year winter
last_nc_yr <- as.numeric(nc_year) - 1
last_yr_anom_path <- str_c(oisst_path, "annual_anomalies/1982to2011_climatology/daily_anoms_", last_nc_yr, ".nc")

# Load previous year as stack to get the full winter
last_yr_anoms <- stack(last_yr_anom_path)

# Join to 2022
anoms_double <- stack(list(last_yr_anoms, anoms_2022))


# Drop December 2022 so its not influencing the Winter of 2021-2022
dec_2022 <- str_c("X", nc_year, ".12")
not_dec21 <- which(str_sub(names(anoms_double), 1, 8) != dec_2022)
anoms_nodec <- anoms_double[[not_dec21]]


# Set up list of year and month combos to use map() for each season
month_key <- list("Winter" = c(str_c(last_nc_yr, ".12"),
                               str_c(nc_year,  ".01"), 
                               str_c(nc_year, ".02")),
                  "Spring" = str_c(nc_year, c("03", "04", "05"), sep = "."),
                  "Summer" = str_c(nc_year, c("06", "07", "08"), sep = "."),
                  "Fall"   = str_c(nc_year, c("09", "10", "11"), sep = "."))
Code
#  Get mean anoms across seasons
season_stacks <- map(month_key, function(mon){
  
  # Get names from the stack - yyyy.mm
  stack_names <- names(anoms_nodec)
  stack_months <- str_sub(stack_names, 2,8)
  
  # layers with correct months:
  in_season <- which(stack_months %in% mon)
  
  # If there is no dates yet just multiply by zero so we can plot a blank
  if(length(in_season) == 0){
    season_mean <- anoms_nodec[[1]] * 0
  } else {
     # Get mean across those months
     season_mean <- calc(anoms_nodec[[in_season]], mean, na.rm = T)
  }
  
 
  # season_mean <- st_as_stars(rotate(season_mean))
  return(season_mean)
  
  
})
Code
# Plotting the seasons
seas_anom_maps <- function(season, lab_years, temp_lim = 5){
  
 
  # Center the color scale with Dynamic Limits
  temp_limits <- c(-1,1) * temp_lim
  temp_breaks <- c(temp_limits[1], temp_limits[1]/2,  0, temp_limits[2]/2, temp_limits[2])
  temp_labels <- str_c(c(str_c("< ", temp_limits[1]), temp_limits[1]/2, 0, temp_limits[2]/2, str_c("> ", temp_limits[2])), "\u00b0C")
  
  # Make map
  seas_map <- ggplot() +
    geom_stars(data = st_as_stars(rotate(season_stacks[[season]]))) +
    geom_sf(data = world, fill = "gray30", color = "white", size = .25) +
    scale_fill_distiller(palette = "RdBu", 
                     na.value = "transparent", 
                     limit = temp_limits, 
                     oob = scales::squish,
                     breaks = temp_breaks, 
                     labels = temp_labels) +
    guides("fill" = guide_colorbar(title = sst_lab, 
                                   title.position = "right", 
                                   title.hjust = 0.5,
                                   barheight = unit(3, "inches"), 
                                   frame.colour = "black", 
                                   ticks.colour = "black")) +  
    coord_sf(expand = F) +
   map_theme() +
    theme(legend.position = "right", 
          legend.title = element_text(angle = 90)) +
    labs(title = str_c(season, " ", lab_years))
  return(seas_map)

  
}
Code
seas_anom_maps("Winter", str_c(last_nc_yr, "-", nc_year, " through: ", Sys.Date()), temp_lim = 2)

Spring

Code
seas_anom_maps("Spring", str_c(nc_year, " through: ", Sys.Date()), temp_lim = 2)

Summer

Code
seas_anom_maps("Summer", str_c(nc_year, " through: ", Sys.Date()), temp_lim = 2)

Fall

Code
seas_anom_maps("Fall", str_c(nc_year, " through: ", Sys.Date()), temp_lim = 2)

2: Regional Temperature Maps

The regional views of 2022’s sea surface temperature anomalies have been loaded and displayed below to visualize how localized differences changed throughout the year.

Code
# Load the bounding box for Andy's GOM to show they align
poly_path     <- region_paths[[params$region]][["shape_path"]]
region_extent <- st_read(poly_path, quiet = TRUE)

# Pull extents for the region to set crop extent
crop_x <- st_bbox(region_extent)[c(1,3)]
crop_y <- st_bbox(region_extent)[c(2,4)]

# # Expand the area out to see the larger patterns
crop_x <- crop_x + c(-2.5, 3.5)
crop_y <- crop_y + c(-1.5, 1)

# Make a new extent for cropping
region_extent_expanded <- st_sfc(st_polygon(list(
  rbind(c(crop_x[[1]], crop_y[[1]]),  
        c(crop_x[[1]], crop_y[[2]]), 
        c(crop_x[[2]], crop_y[[2]]), 
        c(crop_x[[2]], crop_y[[1]]), 
        c(crop_x[[1]], crop_y[[1]])))))

# convert to sf
region_extent_expanded <- st_as_sf(region_extent_expanded)

Annual Average

Code
# Mask the annual average
reg_ann_anom <- mask_nc(ras_obj = ann_anom_ras, mask_shape = region_extent_expanded)

# Color limit for palettes
temp_limits <- c(-5, 5)
temp_breaks <- c(temp_limits[1], temp_limits[1]/2,  0, temp_limits[2]/2, temp_limits[2])
temp_labels <- str_c(c(str_c("< ", temp_limits[1]), temp_limits[1]/2, 0, temp_limits[2]/2, str_c("> ", temp_limits[2])), "\u00b0C")

# Build Map
ggplot() +
  geom_stars(data = st_as_stars(reg_ann_anom)) +
  geom_sf(data = new_england, fill = "gray30", color = "white", size = .25) +
  geom_sf(data = canada, fill = "gray30", color = "white", size = .25) +
  geom_sf(data = greenland, fill = "gray30", color = "white", size = .25) +
  scale_fill_distiller(palette = "RdBu", 
                     na.value = "transparent", 
                     limit = temp_limits, 
                     oob = scales::squish,
                     breaks = temp_breaks, 
                     labels = temp_labels) +
    guides("fill" = guide_colorbar(title = sst_lab, 
                                   title.position = "right", 
                                   title.hjust = 0.5,
                                   barheight = unit(3, "inches"), 
                                   frame.colour = "black", 
                                   ticks.colour = "black")) +  
  map_theme() +
  theme(legend.position = "right", 
        legend.title = element_text(angle = 90)) +
  coord_sf(xlim = crop_x, 
           ylim = crop_y, expand = F) +
  labs(title = str_c("Regional Temperature Anomalies - ", nc_year),
       subtitle = str_c(nc_year, " through: ", Sys.Date()))

Winter

Code
# Mask the Seasons
seasons_masked <- map(season_stacks, mask_nc, region_extent_expanded)

# Plot the seasons
plot_masked_season <- function(season, year_lab, temp_lim = 5){
  
  # Grab Season
  reg_seas_anom <- seasons_masked[[season]]
  
  # Get temp limits
  # temp_limits <- max(abs(values(reg_seas_anom)), na.rm = T) * c(-1,1) # Dynamic Limits
  # Color limit for palettes
  temp_limits <- c(-1, 1) * temp_lim
  temp_breaks <- c(temp_limits[1], temp_limits[1]/2,  0, temp_limits[2]/2, temp_limits[2])
  temp_labels <- str_c(c(str_c("< ", temp_limits[1]), temp_limits[1]/2, 0, temp_limits[2]/2, str_c("> ", temp_limits[2])), "\u00b0C")
  
  
  # Build Map
  seas_map <- ggplot() +
    geom_stars(data = st_as_stars(reg_seas_anom)) +
    geom_sf(data = new_england, fill = "gray30", color = "white", size = .25) +
    geom_sf(data = canada, fill = "gray30", color = "white", size = .25) +
    geom_sf(data = greenland, fill = "gray30", color = "white", size = .25) +
    scale_fill_distiller(palette = "RdBu", 
                     na.value = "transparent", 
                     limit = temp_limits, 
                     oob = scales::squish,
                     breaks = temp_breaks, 
                     labels = temp_labels) +
    guides("fill" = guide_colorbar(title = sst_lab, 
                                   title.position = "right", 
                                   title.hjust = 0.5,
                                   barheight = unit(3, "inches"), 
                                   frame.colour = "black", 
                                   ticks.colour = "black")) +  
  map_theme() +
  theme(legend.position = "right", 
        legend.title = element_text(angle = 90)) +
    coord_sf(xlim = crop_x, 
             ylim = crop_y, expand = F) +
    labs(title = str_c(season, " - ", year_lab))
  
  return(seas_map)
  
}
Code
plot_masked_season("Winter", year_lab = str_c(last_nc_yr, "-", nc_year, " through: ", Sys.Date()), temp_lim = 3.5)

Spring

Code
plot_masked_season("Spring", year_lab = str_c(nc_year, " through: ", Sys.Date()), temp_lim = 3.5)

Summer

Code
plot_masked_season("Summer", year_lab = str_c(nc_year, " through: ", Sys.Date()), temp_lim = 3.5)

Fall

Code
plot_masked_season("Fall", year_lab = str_c(nc_year, " through: ", Sys.Date()), temp_lim = 3.5)

Heatwave Progression

Looking specifically at the last heatwave event, we can step through how the event progressed over time, and developing pockets or warmer/colder water masses.

Code
# Identify the last heatwave event that happened
last_event <- max(region_hw$mhw_event_no, na.rm = T)

# Pull the dates of the most recent heatwave
last_event_dates <- region_hw %>% 
  filter(mhw_event_no == last_event) %>% 
  pull(time)


# Buffer the dates, start 7 days before
event_start <- (min(last_event_dates) - 7)
event_stop  <- max(last_event_dates)
date_seq <- seq.Date(from = event_start,
                     to   = event_stop,
                     by   = 1)



# Load the heatwave dates
data_window <- data.frame(
  time = c(min(date_seq) , max(date_seq) ),
  lon  = crop_x,
  lat  = crop_y)


# Pull data off box
hw_stack <- oisst_window_load(data_window = data_window, 
                              anomalies = T, mac_os = "mojave")


#drop any empty years that bug in
hw_stack <- hw_stack[map(hw_stack, class) != "character"]


##### Format the layers and loop through the maps  ####

# Grab only current year, format dates
this_yr   <- stack(hw_stack)
day_count <- length(names(this_yr))
day_labs  <- str_replace_all(names(this_yr), "[.]","-")
day_labs  <- str_replace_all(day_labs, "X", "")
day_count <- c(1:day_count) %>% setNames(day_labs)

# Progress through daily timeline to indicate heatwave status and severity
hw_timeline <- region_hw %>% 
  filter(time %in% as.Date(day_labs))
Code
####  Plot Settings:

# Set palette limits to center it on 0 with scale_fill_distiller
limit <- c(max(values(this_yr), na.rm = T) * -1, 
           max(values(this_yr), na.rm = T) )


# Plot Heatwave 1 day at a time as a GIF
day_plots <- imap(day_count, function(date_index, date_label) {
  
  # grab dates
  heatwaves_st  <- st_as_stars(this_yr[[date_index]])
  
  #### 1. Map the Anomalies in Space
  day_plot <- ggplot() +
    geom_stars(data = heatwaves_st) +
    geom_sf(data = new_england, fill = "gray30", color = "white", size = .25) +
    geom_sf(data = canada, fill = "gray30", color = "white", size = .25) +
    geom_sf(data = greenland, fill = "gray30", color = "white", size = .25) +
    geom_sf(data = region_extent, 
            color = gmri_cols("gmri blue"), 
            linetype = 2, size = 1,
            fill = "transparent") +
    scale_fill_distiller(palette = "RdYlBu", 
                         na.value = "transparent", 
                         limit = limit) +
    map_theme(legend.position = "bottom") +
    coord_sf(xlim = crop_x, 
             ylim = crop_y, 
             expand = T) +
    guides("fill" = guide_colorbar(
      title = "Sea Surface Temperature Anomaly \u00b0C",
      title.position = "top",
      title.hjust = 0.5,
      barwidth = unit(3, "in"),
      frame.colour = "black",
      ticks.colour = "black")) 
  
  # Set colors by name
  color_vals <- c(
    "Sea Surface Temperature" = "royalblue",
    "Heatwave Event"          = "darkred",
    "MHW Threshold"           = "coral3",
    "Daily Climatology"        = "gray30")
  
  
  
  #### 2.  Plot the day and the overall anomaly to track dates
  date_timeline <- ggplot(data = hw_timeline, aes(x = time)) +
    geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
    geom_line(aes(y = hwe, color = "Heatwave Event")) +
    geom_line(aes(y = cse, color = "Cold Spell Event")) +
    geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
    geom_line(aes(y = mcs_thresh, color = "MCS Threshold"), lty = 3, size = .5) +
    geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = .5) +
    scale_color_manual(values = color_vals) +
    
    # Animated Point /  line
    geom_point(
      data = filter(hw_timeline, time == as.Date(date_label)),
      aes(time, sst, shape = factor(mhw_event)), 
      color = gmri_cols("gmri blue"), 
      size = 3, show.legend = FALSE) + 
    geom_vline(data = filter(hw_timeline, time == as.Date(date_label)),
               aes(xintercept = time), 
               color = "gray50", 
               size = 0.5,
               linetype = 3,
               alpha = 0.8) + 
    labs(x = "", 
         y = "",
         color = "",
         subtitle = "Regional Temperature \u00b0C",
         shape = "Heatwave Event") +
    theme(legend.position = "bottom")
  
  
  ####  3. Assemble plot(s)
  p_layout <- c(
    area(t = 1, l = 1, b = 2, r = 8),
    area(t = 3, l = 1, b = 8, r = 8))
  
  # plot_agg <- (date_timeline / day_plot) + plot_layout(heights = c(1, 3))
  plot_agg <- date_timeline + day_plot + plot_layout(design = p_layout)
  
  
  return(plot_agg )
  
  
})


walk(day_plots, print)

3: Ranking the Rate of Warming

If we look at the rates of change from 1982-2021 for each grid cell, rather than the observed temperature, it is possible to rank how hot each location on earth is warming relative to all the others.

Once we have the rankings, we can then take the average ranking within the Gulf of Maine we can obtain the average warming rank for the area compared to the rest of the globe.

Code
# 1. Warming Rates and Rankings
rates_path <- paste0(oisst_path, "warming_rates/annual_warming_rates")
rates_stack_all <- stack(str_c(rates_path, "1982to2021.nc"), 
                         varname = "annual_warming_rate")
ranks_stack_all <- stack(str_c(rates_path, "1982to2021.nc"), 
                         varname = "rate_percentile")
Code
# Get the rank information that go with the original extent used by the timelines
ranks_masked <- mask_nc(ranks_stack_all, region_extent)
rates_masked <- mask_nc(rates_stack_all, region_extent)
region_ranks <- get_masked_vals(ranks_masked, rates_masked)

# Prep it for text input.
avg_rank <- region_ranks$`Mean Rank` *100
avg_rate <- region_ranks$`Mean Rate`
low_rank <- region_ranks$`Min Rank` *100
low_rate <- region_ranks$`Min Rate`
top_rank <- region_ranks$`Max Rank` *100
top_rank <- ifelse(top_rank == 100, "greater than or equal to 99.5", top_rank)
top_rate <- region_ranks$`Max Rate`

Based on data from 1982-2021, the warming rates of Gulf of Maine have been some of the highest in the world. The area as a whole has been increasing at a rate of 0.044\(^{\circ}C/year\) which is faster than 95.9% of the world’s oceans.

Over that same period locations within the Gulf of Maine have been warming at rates as low as 0.017\(^{\circ}C/year\) and as rapidly as 0.094\(^{\circ}C/year\), corresponding to ranks as low as 60.3% and as high as greater than or equal to 99.5%.

Mapped below are the corresponding warming rates and their global rankings.

Code
# For the full map  we mask again, but zoom out a little

# Mask again using the expanded mask so we can zoom out
ranks_masked <- mask_nc(ranks_stack_all, region_extent_expanded)
rates_masked <- mask_nc(rates_stack_all, region_extent_expanded)


# Make stars object
rank_stars <- st_as_stars(ranks_masked)
rates_stars <- st_as_stars(rates_masked)

Warming Rate Map

Code
####  Rates Map  ####

# Make Contours if desired
# rates_contour <-  st_contour(x = rates_stars, na.rm = T, breaks = seq(0.85, 1, by = 0.02))

# # rates map Original
# rates_map <- ggplot() +
#   geom_stars(data = rates_stars) +
#   geom_sf(data = new_england, fill = "gray90", size = .25) +
#   geom_sf(data = canada, fill = "gray90", size = .25) +
#   geom_sf(data = greenland, fill = "gray90", size = .25) +
#   geom_sf(data = region_extent, 
#           color = "black", 
#           fill = "transparent", linetype = 2, size = 0.5) +
#   scale_fill_viridis_c(option = "plasma", na.value = "transparent") +
#   map_theme() +
#   coord_sf(xlim = crop_x, 
#            ylim = crop_y, expand = F) +
#   guides("fill" = guide_colorbar(
#     title = "Annual Temperature Change \u00b0C/year",
#     title.position = "top", 
#     title.hjust = 0.5,
#     barwidth = unit(2.5, "in"),
#     frame.colour = "black",
#     ticks.colour = "black"))



# Look at the spread
# hist(rates_stars$Annual.Sea.Surface.Temperature.Warming.Rate)

# Reclassify to discrete bins
rates_reclass <- rates_stars %>% 
  mutate(
    warm_rate = Annual.Sea.Surface.Temperature.Warming.Rate,
    rates_reclass = case_when(
      warm_rate > .09  ~ "> 0.1", 
      warm_rate > .08 ~ "0.08 - 0.1", 
      warm_rate > .06 ~ "0.06 - 0.08", 
      warm_rate > .04 ~ "0.05 - 0.06", 
      warm_rate > .04 ~ "0.04 - 0.05", 
      warm_rate > .02 ~ "0.02 - 0.04", 
      warm_rate < 0.2 ~ "< 0.02",
      TRUE ~ "NA"),
    rates_reclass = factor(rates_reclass, levels = c(
      "> 0.1", "0.08 - 0.1", 
      "0.06 - 0.08", "0.05 - 0.06", "0.04 - 0.05", 
      "0.02 - 0.04", "< 0.02",
      "Outside Area Measured"))
  ) %>% select(rates_reclass)


# Testing discrete scales
rates_map <- ggplot() +
  geom_stars(data = rates_reclass) +
  geom_sf(data = new_england, fill = "gray90", size = .25) +
  geom_sf(data = canada, fill = "gray90", size = .25) +
  geom_sf(data = greenland, fill = "gray90", size = .25) +
  geom_sf(data = region_extent, 
          color = "black", 
          fill = "transparent", linetype = 2, size = 0.5) +
  scale_fill_brewer(palette = "YlOrRd", 
                    na.value = "transparent", 
                    na.translate = F,
                    direction = -1) +
  coord_sf(xlim = crop_x, 
           ylim = crop_y, expand = F) +
  map_theme(legend.position = "bottom") +
  guides("fill" = guide_legend(
    title = "Annual Temperature Change \u00b0C/year",
    title.position = "top", 
    title.hjust = 0.5,
    label.position = "bottom", 
    keywidth = unit(1.5, "cm"),
    reverse = T,
    nrow = 1))


rates_map
grid::grid.raster(lab_logo, x = 0.05, y = 0.03, just = c('left', 'bottom'), width = unit(0.8, 'inches'))

Warming Rank Map

Code
####  Warming Ranks Map  ####

# Make Contours if desired
ranks_contour <-  st_contour(x = rank_stars, na.rm = T, breaks = seq(0.85, 1, by = 0.02))


# # Original Map
# # ranks map
# ranks_map <- ggplot() +
#   geom_stars(data = rank_stars) +
#   geom_sf(data = ranks_contour, fill = "transparent", color = "gray30", size = 0.1) +
#   geom_sf(data = new_england, fill = "gray90", size = .25) +
#   geom_sf(data = canada, fill = "gray90", size = .25) +
#   geom_sf(data = greenland, fill = "gray90", size = .25) +
#   geom_sf(data = region_extent, 
#           color = "black", 
#           fill = "transparent", linetype = 2, size = 0.5) +
#   scale_fill_viridis_c(option = "plasma",
#                        na.value = "transparent",
#                        limit = c(0.85, 1),
#                        oob = scales::oob_squish) +
#  map_theme() +
#   coord_sf(xlim = crop_x, 
#            ylim = crop_y, expand = F) +
#   guides("fill" = guide_colorbar(
#     title = "Global Percentile of Warming Rates",
#     title.position = "top", 
#     title.hjust = 0.5,
#     barwidth = unit(2.5, "in"),
#     frame.colour = "black",
#     ticks.colour = "black")) +
#   labs(caption = "Ranking color scale truncated to display ranges of 0.85-1 
#                   Lower values will display as 0.85 or 85%
#                   Contour values every 2%")


# Look at the spread
# hist(rank_stars$Annual.Sea.Surface.Temperature.Warming.Rate.Rank)

# reclassify to discrete bins
ranks_reclass <- rank_stars %>% 
  mutate(
    warm_rank = Annual.Sea.Surface.Temperature.Warming.Rate.Rank,
    ranks_class = case_when(
      warm_rank > .999  ~ "> 99.9%", 
      warm_rank > .99  ~ "99 - 99.9%", 
      warm_rank > .98 ~ "98 - 99%", 
      warm_rank > .95 ~ "95 - 98%", 
      warm_rank > .90 ~ "90 - 95%", 
      warm_rank > .80 ~ "80 - 90%", 
      warm_rank < .80 ~ "< 80%",
      TRUE ~ "NA"),
    ranks_class = factor(ranks_class, levels = c(
      "> 99.9%", "99 - 99.9%", "98 - 99%", 
      "95 - 98%", "90 - 95%", 
      "80 - 90%", "< 80%"))
  ) %>% select(ranks_class)


ranks_map <- ggplot() +
  geom_stars(data = ranks_reclass) +
  geom_sf(data = new_england, fill = "gray90", size = .25) +
  geom_sf(data = canada, fill = "gray90", size = .25) +
  geom_sf(data = greenland, fill = "gray90", size = .25) +
  geom_sf(data = region_extent, 
          color = "black", 
          fill = "transparent", linetype = 2, size = 0.5) +
 scale_fill_brewer(palette = "YlOrRd", 
                    na.value = "transparent", 
                    na.translate = F,
                    direction = -1) +
  map_theme(legend.position = "bottom") +
  coord_sf(xlim = crop_x, 
           ylim = crop_y, expand = F) +
  guides("fill" = guide_legend(
    title = "Global Percentile of Warming Rates",
    title.position = "top", 
    title.hjust = 0.5,
    label.position = "bottom", 
    keywidth = unit(1.5, "cm"),
    reverse = T,
    nrow = 1)) +
  labs(caption = "Pixels ranked globally based on annual warming rate.")

# plot both
ranks_map
grid::grid.raster(lab_logo, x = 0.05, y = 0.03, just = c('left', 'bottom'), width = unit(0.8, 'inches'))

 

A work by Adam A. Kemberling

Akemberling@gmri.org